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For an equation of state in which pressure is a function only of density, the analysis of Newtonian 
stellar structure is simple in principle if the system is axisymmetric, or consists of a corotating 
binary. It is then required only to solve two equations: one stating that the "injection energy," 
K, a potential, is constant throughout the stellar fluid, and the other being the integral over the 
stellar fluid to give the gravitational potential. An iterative solution of these equations generally 
diverges if k is held fixed, but converges with other choices. We investigate the mathematical 
reason for this convergence/divergence by starting the iteration from an approximation that is 
perturbatively different from the actual solution. A cycle of iteration is then treated as a linear 
"updating" operator, and the properties of the linear operator, especially its spectrum, determine 
the convergence properties. For simplicity, we confine ourselves to spherically symmetric models in 
which we analyze updating operators both in the finite dimensional space corresponding to a finite 
difference representation of the problem, and in the continuum, and we find that the fixed-K operator 
is self-adjoint and generally has an eigenvalue greater than unity; in the particularly important case 
of a polytropic equation of state with index greater than unity, we prove that there must be such 
an eigenvalue. For fixed central density, on the other hand, we find that the updating operator has 
only a single eigenvector, with zero eigenvalue, and is nilpotent in finite dimension, thereby giving 
a convergent solution. 



I. INTRODUCTION 
A. Background 

For a star, or a binary pair of stars, modeled as a barotropic fluid rotating with a known dependence of the rotation 
rate f2 on the radial distance zu from the rotation axis, there are two Newtonian forces per unit mass acting on a fluid 
element. The pressure gradient contributes —p~^'S/p where p is the stellar fluid mass density; gravity contributes — V$, 
where $ is the Newtonian potential. The equivalence of these forces per unit mass to the centripetal acceleration is 

Vp + pS/^ ^ pwfl^ {tu) . (1.1) 

For "cold" stellar fluid, as in a neutron star, thermal transport is unimportant and the pressure can be considered 
to be a function only of the density. The force equation p.ip can then be written as 

V ( + $ - -m^n^{w)) = , (1.2) 



where h, the specific enthalpy, is defined by 
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h{p):=r^. (1.3) 



P 



The equations of stellar structure can then be taken as 



^-^w^n^{w) = K (1.4) 
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p{x') 




(1.5) 
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in which the second equation is the Poisson integral for the gravitational potential With h a specified function 
of p, these two equations constitute the basis for finding ^{x) and p{x), and hence for solving the problem of stellar 
structure. 

Two essentially different iterative methods have been used to solve this nonlinear system; each was developed first 
in the Newtonian context and then extended to full general relativity: A Newton- Raphson method was first used by 
James [ij and Stoeckly [2| (see also [1, 3 ) • In this method the two equations (|1.4p , (|1.5p must be solved simultaneously, 
requiring the inversion of a large matrix. The second method, the so-called self-consistent field (SCF) method, first 
introduced in the context of rotating stars by Ostriker and Mark Q , is much simpler computationally. In this method 
one first solves one of the pair (|1.4p . (|1.5p and then the other. 

In the SCF (as well as the Newton-Raphson) method, iteration requires fixing two parameters to be held constant 
during the iteration. Typically one of these is the rotation law, such as the condition of uniform rotation ~ constant. 
The most obvious choice for a second condition would be a fixed value of the injection energy k in Eq. p.4p . For this 
choice, the iteration starts with a guess for a density profile p{x); Eq. (|1.5p can then be solved for ^{x); this result 
for ^{x) is put in Eq. p.4p which is solved for h{x); finally, from the form of the functional relationship between p 
and h, a new, iterated, solution is found for p{x). This cycle is then repeated. 

It is found that with k held fixed, the SCF iteration does not converge. For other choices, however, the iteration 
does converge. For rotating stars, for example, the SCF iteration converges if the central density is held fixed. To 
achieve convergence, Ostriker and Mark (and many subsequent authors) fix the total mass, while Hachisu jB] fixes the 
ratio of polar to equatorial radius. 

Although computational astrophysicists have found success using the SCF method, no explanation has ever been 
given for the relationship between convergence of iteration and the fixed conditions of the iteration. Ostriker and Marck 
[5| mention unpublished work showing stability for a particularly simple case (the n = 1 polytrope to be discussed 
below) but we are not aware of other analytic studies of iterative stability of the SCF method of constructing stellar 
models. 

In this paper we provide a mathematical explanation for the convergence/divergence properties of the SCF iteration. 
Part of the motivation for doing this is the hope that improved understanding of the SCF method will lead to simpler 
schemes for computing the structure of rapidly rotating neutron stellar models in general relativity. But most of the 
motivation is curiosity about a feature of computations that has been known but unexplained for more than forty 
years. 

We have approached this problem by considering the iteration to start very close to a solution of Eqs. (|1.4p . (jl.Sp . 
That is, we start with an initial guess for p(x) that is perturbatively different from the exact function that would 
solve Eqs. (|1.4p . (jl.Sp . The exact structure problem (|1.4p . (|1.5p is then replaced by the equations for Sh, S^, 6k, 
Sp, the differences between the true values of h, $, k, p and the values at the start of an iteration cycle. Higher order 
terms in these perturbations are dropped, resulting in linear equations for the perturbations 



with the rotation law fixed during the iteration, Jfi = 0. A cycle of iteration then starts with an initial perturbation 
d^, and ends with an iterated perturbation S^. One round of iteration, therefore, constitutes a linear "updating" of 
the density perturbation 5p{x). Whether or not this linearized iteration converges depends on the properties of the 
linear operator that performs the updating. 

It is fairly clear that convergence of the linearized problem of Eq. (|1.6p is a necessary condition for convergence of 
the SCF iterative solution to the exact equations (|1.4p . (|1.5p . Our working assumption is that it is also a sufficient 
condition; i.e. , if the linearized iteration converges, then iteration of the exact problem will also converge if the 
iteration is started close enough to the exact solution. Our numerical results support this assumption: exact iteration 
converges if and only if linearized iteration converged. 

The analysis of convergence then becomes a study of the properties of the linear updating operator. We find that 
for fixed-K iteration the linear updating operator is self-adjoint and (not surprisingly) generally has an eigenvalue 
larger than unity. Less expected is the result of the study of the cases in which iteration does converge. Here it is 
found that the eigenbasis is not complete, but (for a finite difference representation of the equations) is nilpotent, and 
can usefully be understood with a Jordan decomposition into generalized eigenvectors. 

The evidence for these explanations, for both rotating stars and binaries, is extensive, but consists largely of 
numerical results. The simplicity of the spherical case, on the other hand, allows some very definitive and interesting 
mathematical results and physical insights, and it is that case that we present here. 





(1.6) 
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B. Spherically symmetric equations 



In spherical symmetry the structure problem is that of determining $(r) and p(r), where r is the radial distance 
from the center of symmetry. Since 51 = in the spherical case, the equation of hydrodynamic equilibrium reduces to 



/i + <f> = K . 



(1.7) 



It is convenient to consider the enthalpy /i, rather than the density, to be the unknown structure function, so we 
assume that p is a known invertible functional of h given by 



P = Q[h] ■ 

In the particular, and very common case of a polytropic equation of state, 

p = Kp^+^/" , 

the enthalpy function is 

h = i^(l + , 

so that 

g{h) = [h/K{l + n)r . 
The Poisson equation for spherical symmetry, in terms of h, is 



g[h{r')y^ dr' + I Q[h{r')]r' dr' 



(1.8) 
(1.9) 
(1.10) 
(1.11) 

(1.12) 



and we assume that p is finite at r = 0, so that $(0) is finite. With differentiation this system can be cast as the 
following differential equation for h(r): 



1 



AirGr^ dr \ dr 



r'^]^-g[h{r)]. 



(1.13) 



Equation (|1.12p is equivalent to this differential equation with the constraint that h{r) be analytic at r = 0. The 
value r = i? at which h first vanishes determines the radius R of the stellar surface and the boundary of the domain 
in which Eq. (|1.7p applies. 

The linearization of the spherical problem leads to the equations 



6h + 5^ = Sk 



(1.14) 



5$ = -AttG 



V{r')Sh(r'r'^r'^dr' + / V{r')5h{r'r'^ r' dr' 



in which 



(1.15) 



V := dg/dh. 



(1.16) 



In Eq. (|1.15p we have have ignored changes in the stellar radius R, because their contribution can be shown to be of 
order higher than linear, as long as the star is compressiblc"'^ . The basis of our analysis is Eqs. (|1.14p and (|1.15p . with 
different choices of the form of V, and different choices of what is held fixed during iteration. 



For a polytrope of index n, the contribution in Eq. 111.151 1 of the change 5R in radius can be shown to be of order 1 + n in SR. 
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C. Outline 



Our study of the SCF solution of spherical structure was guided by, and aimed at explaining the known numerical 
phenomena in nonspherical SCF iteration, along with our own numerical discoveries. These come almost entirely from 
computations on a finite difference grid with barotropic equations of state and include the following, (i) Iteration 
with fixed k appears always to diverge, (ii) For rotating stellar models, iteration with fixed central density always 
converges. 

For spherical stellar models, by analyzing the mathematical properties of the linearized updating operator (as 
opposed to studying its numerical results), we have been able to show the following: (i) For iteration with k fixed, 
the linearized updating operator is self-adjoint and there is a complete basis of density perturbations. (This set of 
perturbations is not related to the likewise complete and orthogonal set of density perturbations that are eigenmodes 
of radial oscillations of spherical stellar models.) For models with polytropic index n we have shown that there is a 
single eigenvector with eigenvalue greater than n, and that all others must have eigenvalues less than n. This discovery 
motivated a study of models with n < 1, and it was found that, indeed, for sufficiently small (and astrophysically 
irrelevant) n, the fixed-K operator does have all linearized eigenvalues less than unity, and iteration does converge, 
(ii) For fixed central density on a finite difference grid, there is only a single eigenvector and the single eigenvalue zero. 
This spectrum has a very clear physical explanation. For the continuum version of the linearized updating operator for 
fixed central density, we show that iteration must converge, and that there can be no bounded eigenvector. The single 
eigenvector of the discrete implementation corresponds, in the continuum, to a delta function, (iii) For iteration 
with fixed density at some radius r/ > that is less than the stellar radius i?, we show that there is an infinite 
number of eigenvectors, that these eigenvectors do not form a complete basis, and that the updating may or may not 
converge. In the finite difference implementation of this problem, with N grid zones in the interval (0, i?), the number 
of eigenvectors is equal to the number of grid zones in the interval (0,r^). 

The remainder of this paper is organized as follows. In Sec. |TT] we present the specialization of the spherically 
symmetric iteration scheme to the fixed-K case. We show that the linearized updating operator is self adjoint (for 
a suitably chosen function space and inner product). We derive a strict bound on the eigenvalues in the case of a 
polytropic equation of state, and (less useful) bounds for more general equations of state. We turn, in Sec. Illli to 
iteration with density fixed at some radius r/. Finite difference results are given that show for = that there is only 
a single eigenvector, with zero eigenvalue. For < r/ < R, we show that the number of eigenvectors is proportional to 
the choice of r/. We give a physical explanation of these numerical results and then consider the equivalent problem 
in the continuum. We show that for rf — the iteration must converge, and we show that there can be no bounded 
eigenvector. We go on to show that for rf ^ the iteration problem has aspects both of the fixed-K problem and of 
the fixed central density problem. The paper is summarized, and conclusions given, in Sec. IIVI 



II. FIXED-K ITERATION 



A. Self-adjoint linearized updating operator 

If we set Sk^O in Eqs. (|1.14p . (jl.lSp then Sh = -S^ and we have 



^/i"™(r) = 47rG 



(2.1) 



where is the linear updating operator for fixed n. The inverse of this linearized updating operator is found, by 
differentiating Eq. (|2.ip , to be 



L-\v{r))=~ 



1 



, dv 



AnGr^Vir) dr \ dr 



(2.2) 



Any smooth function v{r) that results from an application of in Eq. (|2.ip will have the property at the stellar 
surface r = i? that 



dv 
dr 



r=R- 



r=R^ 



(2.3) 



We take this to be one of the conditions on the function space on which and ^ operate. The second condition 
is that v{r) must be finite at the stellar center r = 0. 
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As our inner product on this space we choose 



from which we get 



Vi ■ V2 ~ r 'P{r)vi{r)v2{r) dr . 
Jo 



(2.4) 



4nG 



vir 



, dv2 
dr 



V2r 



, dvi 
dr 



(2.5) 



With the conditions that the functions vi , V2 are well behaved at r = and satisfy Eq. (|2.3|) , the right hand side 
above vanishes and we conclude that is self-adjoint^. We shall see that L^^ has no zero eigenvalues, and hence is 
invertible, and L„, its inverse, is self adjoint. (We note here that for nonspherical models it is equally simple to prove 
that the fixed-K updating operator is self-adjoint.) 



B. Polytropes and eigenvalue bounds 

We now specialize to the polytropic equations of state. In this case, the g function in Eq. is used in Eq. (|1.13p . 

To simplify notation we follow common convention 7] and introduce dimensionless variables 

ri/2 



in which ho is the unperturbed value of h at the stellar center r = 0. In terms of these variables the nonlinear equation 
of structure Eq. (|1.13p becomes the Lane-Emden equation^ 

For fixed-K perturbations about a solution of Eq. (|1.12p . the eigenvalue problem £«(«) — for Eq. (|2.ip is 
equivalent to solving the inverse problem L~^{v) = (1/A)w for the differential operator in Eq. (|2.2p . With 
V = dg/dh = nh"~^ /[K{1 + n)]", and with the dimensionless radial variable ^ this becomes 

1 d f ,-:>dv\ n f h \ "^^ 

We next introduce the notation / = ^0 for the Lane-Emden equation, and F = for the perturbation equation, 
and we rewrite these equations, respectively, as 



/ (2.9) 
F . (2.10) 



In the first of these, /(^) is to be considered an unknown function to be found on < ^ < ^max by solving the equation 
subject to the normalization / = ^ -|- 0{^^). In the second equation //^ is considered to be a known function of ^, 
and the eigenequation is to be solved subject to the normalization F = + 0{£^^) and 

rfWle=Uax=0' (2-11) 

which is equivalent to Eq. ()2.3p . 



^ We use "self-adjoint" in tlie pliysicist's sense of an operator symmetric with respect to the inner product defined above, when the 
operator is restricted to smooth functions. We do not characterize its domain. 
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With a solution for (//C)""^ = 6"" 
as particular cases of the equation 



G(^) treated as a known function, we can view both Eqs. (|2.9p and (|2.10|) 



2 - 'kGiO^, 



(2.12) 



with normalization (j) — + C(C ) ■ For this equation we know that if fc = 1, then (j) vanishes at ^ = ^max and is 
positive for ^ < ^max- In Fig. [TJ solutions of Eq. (|2.12p are given in which G{£_) corresponds to the particular case 
n = 2. The curve (j)Q shows the solution for fc = 1, i.e., the solution corresponding to the equilibrium profile of the 
star. The figure also shows eigenfunctions, solutions corresponding to condition (|2.1ip . Shown are the eigenfunction 
for the lowest eigenvalue, 02 for the next higher eigenvalue, and 03 for the next. It seems intuitively clear that 
the eigenvalue k for 0i must be less than unity, since Eqs. (|2.1ip and (|2.12p require that 0i be "less curved" than 
00- Similarly, the eigenvalue for 02, and all other eigenfunctions, must be larger than unity. We now prove that this 
must be so. 

We start by proving that the zeroes and extrema of a smooth solution to Eq. (|2.12|) must alternate. Consider two 
extrema of a solution. There must be a point between those two extrema at which cPip/d^'^ = 0. But Eq. (|2.12p 
requires that = at that point. Between any two zeros, of course, there must be an extremum. Thus zeroes and 
extrema alternate, as claimed. We next consider two functions 0a and 0s that are solutions of Eq. (|2.12p . with the 



1.5 - 

— t 




FIG. 1: Various solutions of (j>,(i = — G(^) for n — 2. 

starting condition (f) = ^ + 0{£,^). Let these two solutions (not necessarily eigenfunctions) correspond respectively to 
kA and ks, with fc^ > kB- Suppose that there is some point ^ = a such that 0a > 0s > for < ^ < a. We then 
consider 

{(l)B - (I)a)^^^ = kA(l)A - kB(l)B for < ^ < a . (2.13) 

By our assumptions, the right hand side is everywhere positive. But the function 0^ — (t>A starts with value zero at 
^ = and with a zero derivative. It follows from Eq. (I2.13P that 0_b — 4>a must be positive for ^ < a, which contradicts 
our assumptions that (pA > 0_b ^ for < ^ < a. This proves that there can be no interval < ^ < a on which 
4' A > 05 > 0. Essentially the same argument shows that 0b is positive at the first zero of 0a, and more generally 
that the ^ value of the first zero of a solution of Eq. p.l2p (subject to the boundary conditions at ^ = 0) decreases 
as k increases. This immediately confirms that an eigenfunction, like 0i, with no zero, must correspond to a value of 
k less than unity. We also conclude that 02 and 03, and any eigenfunction that has extrema intermediate between 
and ^max, must have a zero for ^ < ^max, and hence an eigenvalue k that is larger than unity. This completes the 
proof that there is one and only one eigenvalue k that is smaller than unity. 

When applied to the problem of Eqs. (|2.7p and (|2.8p . this tells us that there is one and only one eigenvalue A of 
the updating operator that is larger than the polytropic index n. For astrophysically relevant models, which have 



7 



n > 1, this guarantees that fixed-K iteration will have an eigenvalue of the updating operator that is larger than unity, 
and that updating will not converge. It suggests, but does not guarantee that there will only be a single updating 
eigenvalue A that is larger than unity. This, however, does turn out to be what we have found in numerical studies 
(see below). 

Though a polytropic equation of state with ri < 1 is not astrophysically plausible, it is interesting since the above 
analysis shows that fixed-K iteration need not diverge for n < 1. We have found numerically that for < n < 0.016 
the nonlinear fixed-K iteration is, indeed, convergent. 

A particularly simple example of the above analysis is the case n — 1, for which the Lane-Emden equation (|2.7p 
admits an analytical solution 

e(0 = ^ (2.14) 

which vanishes at the surface ^max — tt. In this case the eigenvalue problem of Eq. (|2.8p is a spherical Bessel differential 
equation 

and is analytically tractable. The solutions regular at the origin are zeroth order spherical Bessel functions: 



v^'HO=a. ^^^^^^ (2.16) 



where k numbers the eigenfunctions and the ak are normalization constants. To satisfy the boundary condition (j2.3p . 
the eigenvalues are required to have the values 

Afe = (fc-i)"', fc = l,2,.... (2.17) 




FIG. 2: For n = 1 the linearized problem has simple closed form solutions to the problem illustrated in Fig.[T] In this case the 
solution of the Lane-Emden equation gives ^max = tv and 00 ~ sin^. With the eigensolutions (f>k scaled to have = 5 + C(5^) 
the first three eigenmodes are 01 = 2 sin (5/2); 02 ~ (2/3) sin (3^/2); and 03 — (2/5) sin (5^/2). 



The resulting functions 0fc(^) :— S.v^'^^S,) are plotted in Fig.O It is simple to check that the eigenfunctions satisfy 
the orthogonality condition 

^W.„(fc')=/' ^\(fc)(^)„(fc')(^)d5 ^0 iffc^fc'. (2.18) 



8 



For a nonpolytropic equation of state, we can arrive at somewhat weaker results for bounds on the eigenvalues. 
From Eq. (|1.13p and the cigenproblem associated with its linearization we have 



2 -Af (2.19) 



Here / = rh, F = rSh and 



dr 

cPF 1 

^ - -jBF. (2^20) 



A = ATrGg/h B -AwGdg/dh. (2.21) 



Equation (|2.19p is solved and r„iax is defined as the first zero of the solution. That solution is then used in B, so that 
Eq. (|2.20p is considered to be a linear cigenproblem for F. The boundary condition on that equation is dF/dr = at 

^ ^max ■ 

With only slight modification, the argument used in the polytropic case can be used to show: (i) If there is a 
constant ci such that B > ciA for all r S (OjTmax), then the largest A must be greater than ci. (ii) If there is a 
constant C2 such that B < C2v4 for all r G {0,r„^ax), then all As except the largest, must be less than C2. 



C. Numerical investigations 

In terms of the polytropic variables of Eq. (|2.6p . the updating equation (|2.ip can be compactly written in the 
dimensionless form, 

5/i--(^) = n / Q{Cr-' e ^ L^iSh"'") , (2.22) 

Jo max(f,,^') 

which motivates the eigenvalue problem 

L.w e((r-^-^^ei(' = Mo. (2.23) 



This integral form is equivalent to inversion of the differential equation (|2.8p under the boundary conditions (|2.3p . 
The stability of the linearized iteration (|2.22p depends on the spectrum of the linear operator L^.. The spectrum is 
obtained by solving Eq. (|2.23p . which requires that the background (unperturbed) solution 8(^) be known. 

Analytical solutions to the Lane-Emdcn equation fi\ exist only for polytropic indexes n = 0, 1,5; for other values 
of n, we have solved Eq. (|2.7p numerically for 6(^), using a predictor-corrector Adams method of adaptive stepsize 
and order. The eigenvalue problem of Eq. (|2.23p is discretized on the grid of radial coordinates ^ of Eq. (|2.6p using 
N equidistant grid points in the star interior: 

= (z-i)A^, i = l,...,N, (2.24) 

where 

AC = ?max/iV (2.25) 

is the grid spacing. Upon discretization, Eq. (|2.23p becomes: 

N 

J2i^-h^j = (2-26) 

where Vi = v{^i) are the components of a vector v formed from the values of the eigenfunction at the grid points, and 

(L«,),, =ne(e,r-^ i-TT^^- (2.27) 

max(^,;,^j) 

The matrix Sij = l/max(^i, ^j) is symmetric, which means that {Lf^)ij has the form 

N 



{L^),, = S^kDkJ (2.28) 



k=l 
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where D is a diagonal matrix with positive entries. Since the diagonal elements of D are positive, the matrix D^^"^, 
defined by D^/^D^/^ = D, is real and has inverse D"-'^/^. The nonsingular similiarity tranformation D^/^L^D"^/^ = 
D^/^(SD)D~^/-^ — D^/^SD^^^ symmetrizes L«;, proving that L„ has a complete basis of real eigenvectors. Iteration 
will converge if and only if all the eigenvalues have magnitude less than unity. 

Once the background solutions 6(^) are known from a numerical solution of Eq. (|2.7p . their values are interpolated 
to the grid points (|2.24p , and the numerical computation of the eigenvalues of the matrix (|2.27p for various values of 
n S (0, 5) is straightforward. (For n > 5 the value of Q never goes to zero; that is, the stellar model it represents 
has infinite radius. In this case, Cmax does not exist, and the eigenvalue problem is not defined. For n = the stellar 
fluid is incompressible and the background solution is not smooth at the stellar surface.) For sufficiently large N, 
these eigenvalues should approach those of the continuum operator. The results of such a computation are plotted in 
Fig. [3l and confirm that one and only one eigenvalue A is greater than n, while all others are lower than n. 

To check a nonpolytropic equation of state we used 

g{h) = {ah + bh^f/^. (2.29) 

which approximates white dwarf equations of state As one might expect, for a ^ bhQ the unperturbed solution 
approaches an n = 3/2 polytrope, while for a <C bho the unperturbed solution approaches an n = 3 polytrope. It is 
also known Q (page 430) that generic solutions for the equation of state (|2.29p . with arbitrary values of bho /a, are 
bounded by the aforementioned polytropes. Upon computation, we found that the eigenspectrum of also exhibits a 
similar behavior. The eigenvalues of this problem vary monotonically between the n = 3/2 and the n = 3 eigenvalues, 
as bho/ a varies from to oo. One may thus be able to infer the nature of the eigenspectrum - and convergence vs. 
divergence - by considering polytropic equations of state that in some sense "bound" the actual equation of state. 




Polytropic index n 



FIG. 3: Eigenvalues A for various polytropic indices n. The maximum eigenvalue Ai is separated from A2, A3, ... by the separatrix 
A = n (dashed line). Black dots denote the eigenvalues (|2.17|l that correspond to the modes of Fig. O 



III. FIXED DENSITY 
A. Fixed central density and finite difference computations 

As a specific case of conditions that lead to convergent iteration, we consider fixed density at some radius. Here 
we start with the simplest case: fixed density at the stellar center. We have studied convergence of the SCF iteration 
for spherically symmetric polytropic models with a wide range of indices. To check a nonpolytropic equation of state 
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we again used Eq. (|2.29p . In all cases we found that the iteration converged. We now analyze why this is so by 
considering small deviations from the solution. 

For spherical models with fixed central density, Sh = —S^ + 5^\r=Q in Eqs. (|1.14|) . (jl.lSp . and we get 



5/i"™(r) = inG [ r{r')5h{rT'' ( ^ - -\ r'^ dr' = Lee„t(<5/^°"^) . (3.1) 

Jo Vmax(r,r') r' J 

Although the updating operator Lccnt is superficially similar to in Eq. (|2.1|) . the operator Lcont is not self-adjoint 
and, as we shall now demonstrate, its eigenspectrum is dramatically different from that of L^. 

The properties of this updating operator become particularly transparent upon discretization. With a discretization 
that is a modification of that in Sec. Ill G\ 

Ar = R/N, -(j-i)Ar, j = l,...,^ (3.2) 

the eigenvalue problem for Lccnt becomes 

N 

^{Lcent)ijVj = AWi (3.3) 
3 = 1 

with 

(Lccnt)^, = 4nGV{r,) ( ^ ^ - -) rf Ar . (3.4) 

Vmax(ri,rj) rj J ^ 

We note that (LcGnt)y = for ^ Vj or, equivalently, i ^ j. This means that the matrix Lcont G M.^^^ is strictly 
lower triangular, that is, a lower triangular matrix with zeroes on the diagonal. It follows that 

det(Lccnt-AI) = A^ = 0, (3.5) 

so that the only eigenvalue can be zero. Below the diagonal, no element is zero, and each column has a different 
length of nonzero entries. (Note that neither r = nor r = i? is included in the grid, so Virj) is always nonzero.) 
It follows that there are A'^ — 1 linearly independent columns, and hence that the rank of the matrix is iV — 1, and 
therefore that there is only a single zero eigenvector. It is easily seen that, modulo scaling, this eigenvector is 



3 

since this satisfies Eq. (|3.3p . with vanishing eigenvalue: 



v\^^ = 6jN, (3.6) 



N N 



(L cent) ijVj^^ = ^(Lccnt)ij'5jAr — (Lccnt)iJV = 0. (3.7) 



The convergence of the linearized iteration is obvious from the strictly lower triangular nature of {Lccnt)ij- When 
this matrix is applied to any column vector, the result is a column with leading entry zero. A second application 
gives zero for the first two elements of the column, etc. It is obvious that the matrix is therefore nilpotent of index 
N. Not only is iteration with (Lcont)y convergent, it reduces any initial perturbation to zero after N iterations. This 
suggests why the SCF method of solving may be so successful. 

Though a complete eigenbasis does not exist, one can construct a basis of generalized eigenvectors, by putting 
the matrix for Lccnt in Jordan canonical form [1, an approach that will be useful for nonspherical models. In 
this block-diagonal form, the subspace corresponding to each block contains a single eigenvector. In our spherically 
symmetric case there is only a single eigenvector in the whole space, so the Jordan canonical form consists of a single 
block. The basis vectors in the canonical form, v'^'^-' with k ~ 1 . . . N , are the Jordan generalized eigenvectors, and 
satisfy 

(Lccnt - AI)v('=) = Lccntv(^) = V^^-D (3.8) 

for k = 2 . . .N, along with the equation for the true eigenvector LccntV*-^-* = 0. The matrices for Lccnt and v^*^^ in 
this basis have the forms 

(Lccnt),, = <5,j-l (^e^)), (3.9) 
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From Eq. (|3.8p it is clear that the apphcation N times of Lccnt gives zero for any of the basis vectors, and hence for 
any vector. This again shows that the operator is nilpotent of index N. 

The spectrum of Lccnt admits a beautifully simple physical interpretation. The first generalized eigenvector, the 
true eigenvector, corresponds to 6h ^ only in the outermost shell. In an iteration cycle, we first solve for the 
potential inside this outermost shell and find that the only change is that the potential is uniformly changed by a 
constant, (5$. Since the central density, and hence the central enthalpy, is kept fixed, we adjust k in Eq. (|1.14|) so 
that Sh = at the origin. But ^$ is the same everywhere in the stellar interior, so by setting Sh to zero at the origin, 
we set it to zero everywhere. Thus an initial perturbation only in the outermost shell is made to vanish in one cycle 
of iteration. This is the physical picture of the mathematical fact Lccntv^^-* = 0. 

The second, generalized eigenvector (v^^^ in the notation of Eq. (|3.8p ). consists of only the two outermost grid zones 
having 5h ^ 0. This means that all zones interior to the outermost zone will have the same change (5$, while the 
outermost zone will have a different value of 5$. By a minor variation of the previous argument we can see that one 
cycle of iteration will eliminate the density perturbation except in the outermost shell, in other words, will convert 
v^^^ to v*^^^. The extension of this viewpoint explains all the generalized eigenvectors, with v^'^^ having Sh ^ only 
in the outermost three zones, and so forth. (Note: For all the generalized eigenvectors except the true eigenvector, 
we could set the density in the outer shell to zero; since density only in that shell is the zero eigenvector, changing it 
doesn't affect the action of the updating operator.) 



B. Fixed central density and the continuum 



The motivation for this paper is primarily the convergence of numerical methods, so the considerations above for 
finite difference computations suffice in practice for fixed central density iteration. As a matter of principle, however, 
it is interesting to consider the continuum equivalent of the finite difference problem of the previous subsection. 

We start by rewriting Eq. (j3.ip in a notation for finding the p + 1 iterant from the pth. 



V{r') 



max(r, r') r' 



Sh{r')^P^ dr' = -AttG J r{r')r' (^1 - ^ j 5h{r')^P'^ dr' . 

(3.10) 

The factor in square brackets in the integral is nonnegative. We let ^ftSax be the maximum on the interval (0, R), of 
the initial deviation |(5/i(r')^'^^ | from the solution and T'max the maximum of ^(r) on (0, i?). We then have 



This inequality and Eq. (|3.10p then gives us 

\6h('^\r)\ |47rGP„ 

and 



(2p+l)! 



,2 IP 



(2p+l)! 



(3.11) 



(3.12) 



(3.13) 



But G'P /{2p + 1)! ^ as p ^ cxo for any finite C, hence the iteration defined by Eq. p.lOp converges. Unlike the 
discrete case, the continuum operator is not nilpotent, since the results of operating on a function a finite number 
of times does not give zero. The name "quasinilpotent," however, is sometimes applied to an operator like icont for 
which the spectrum consists only of zero. 

We can also inquire about the eigenvector problem for the continuum 



47rG 



L 



r{r'y 1 



v{r') dr' = Lccnt(w) = Xv . 



(3.14) 



If we assume that v{r) is bounded, and A ^ then the same argument used to arrive at Eq. (|3.13p tells us that when 
Lccnt is applied n times to v we get 



\v\< 



47rGPn.axi?' 



(2p+l)!' 



(3.15) 
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which vanishes as p —>■ oo showing that no bounded eigenfunction with A ^ can exist. 

We next show that no bounded eigenfunction with A = can exist. Since 'P(r')r'(l — r'/r) is nonnegative in the 
integrand of Eq. (|3.14p . v{r) must change sign in the integral. Let us assume that, after r — 0, the eigenfunction v{r) 
has its smahest zero at r = a. For definitiveness we take v to be positive on (0, a). From Eq. p.l4p . we have 



Xv{a) = -AttG 



1 - - V[r')v{r')dr' 
a ' 



(3.16) 



In the integrand the factor V{r'){l — r' /a) is nonnegative and not identicaUy zero, and by hypothesis v is nonnegative 
and not identically zero. It follows that via) cannot be zero. Since this contradicts our assumption about a zero at 
a, we conclude that v{r) cannot have a zero in (0, i?), and hence a bounded eigenfunction cannot exist. 

If we relax the condition that the eigenfunction must be bounded, and allow distributional solutions, we immediately 
see that the delta function v = 5{r — R) is an eigensolution since it satisfies 



icc„t(<5(r-i?)) =47rG 



' _ _ 1 V{r')5{r' -R)dr' 



= 0. 



(3.17) 



This eigensolution, of course, is the continuum equivalent of the "outer shell only" eigenvector of the discrete problem. 

It is interesting to consider the analog in the continuum of the Jordan canonical form[l3|. This would require a 
definition of the functions that constitute our Banach space on which the operator £cont operates. Without going 
into such detail we can make some interesting observations about a Jordan-like decomposition for £cent- To start we 
note that in an N dimensional context the Jordan basis (for a single zero eigenvalue Jordan block) can be constructed 
starting with the eigenvector v^^^ and proceeding with an inverse operator. (This inverse, of course, is not 

unique, but we can choose it always to give a result orthogonal to v^^\) With this inverse we construct 



,(2) 



,(3) 



^cent 



{V 



If we attempt to follow this pattern in the continuum we can use the inverse of Lcont to be 



^cent(/) 



1 



47rGP(r)r2 dr \ dr' 
With Tj*^^^ = 6{r — R), the analogous sequence of distributions is given by 



(3.18) 



(3.19) 



(3.20) 



Such a sequence - technical objections aside - would give a basis with the property in Eq. p.Sp . The technical 
objection, of course, is that the eigenfunction is a delta function, so that our sequence would consist of more and more 
singular generalized functions. Worse, this sequence in no way resembles the finite dimensional Jordan basis in which 
each subsequent basis vector is an outer shell that is "thicker" than the previous basis vector. 
A more interesting sequence consists of the functions 



,(0) 



1 



,(-1) 



.(-2) 



(3.21) 



This sequence formally satisfies the Jordan basis criterion in Eq. (|3.8p for k = 0,-1, —2, • • • . The limit of this sequence, 
"^(-oo)>5 should in some sense represent the single eigenfunction. That is, w*^^'') should approach the delta function 
at the stellar surface as /c ^ oo. 

We let a simple example suffice to show that, in a rough sense, this is the case. For the n — I polytropic equation 
of state we have from Eqs. (|l.lip and (|1.16p that P = 1/2K, so that the procedure of Eq. (|3.2ip gives 



.,(0) 



.,(-1) - 



-27rG 
K 



r 
31 



-27rGy 



,,(-p) 



K J 5\ 



Y 



r2p 



-27rG 



K J {2p + iy. 



(3.22) 



In intuitive accord with the finite dimensional case, and with the physical picture, as p — > cxd the generalized eigen- 
function, in a rough sense, approaches a density profile that is concentrated at the outer boundary. 
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C. Fixed intermediate density and finite difference computations 



We now consider the case of spherical models with density fixed at distance rf ^ R from the center. With 6h\rj- = 
in Eqs. (|1.14p . (|1.15p we arrive at 



rir') 



J\old 



max(r, r') 



r'2 dr' 



Vir ) —r dr 

max(r/, r) 



Lrf {5h' 



old\ 



(3.23) 



We discretize as in Eq. (|3.2[) . and for convenience we choose rj — (/ — 5) Ar where / is a positive integer < N. In 
place of Eq. ()3.4p we now have 



We notice that this matrix has the structure 



1 




1 




niax(r/, Vj) 


SD 









T 





r?Ar. 



(3.24) 



(3.25) 



where: (i) SD is a (/ — 1) x (/ — 1) square matrix consisting of a symmetric matrix right multiplied by a diagonal 
matrix; (ii) T is a strictly lower triangular {N — / + 1) x {N — / + 1) square matrix; (iii) M is a {N — / + 1) x (/ — 1) 
matrix. To discuss eigensolutions we write column eigenvectors in the form 



u 
w 



where u is a column of length / — 1 and w is a column of length iV — / + 1 . so that 



SD 







u 




SD 


u 


M 


T 




w 




M u + 


T w 



For vectors with u = the eigenproblem reduces to 

T w = Aw . 



(3.26) 



(3.27) 



(3.28) 



Since T is strictly lower triangular we have, from the discussion in Sec. IIII A( that the only eigenvalue is zero, and 
that there is only a single eigenvector. The rest of the — / + 1 dimensional space on which T operates is spanned 
by generalized eigenvectors, as in Sec. IIII Al 

We next consider solutions of the (/ — 1) x (/ — 1) problem 



SD u = Au. 



(3.29) 



We have seen in Sec. Ill Ci that the eigenvectors for this problem are complete in the (/ — 1) x (/ — 1) subspace, hence 
there exist / — 1 eigenvectors in the (/ — 1) x (/ — 1) sector with, in general, distinct eigenvalues. Let u'^'^^ with 
fc=l,2,---/ — 1 represent the set of these column vectors of length / — 1, and let Afc represent the corresponding 
eigenvalues. 

The next step is to define w'^'"') to be the solution of the (A^ — / + 1) x (A^ — / + 1) matrix equation 



(T- Afcl)w 



(k) 



-MwC^) fc = l,2---/-l 



(3.30) 



with I the unit matrix. If we assume that SD has no zero eigenvectors (which is true for all models numerically 
checked) then T — A^I is invertible, since the only eigenvalue of T is zero. This guarantees that solutions of Eq. (|3.30p 
exist for all (A^ — / + 1) values of k. The column vectors combining u^'^^ and w*^*^) are then (A^ — / + 1) eigenvectors, 
since 



(3.31) 



These / — 1 eigenvectors, along with the single zero-eigenvalue eigenvector, are the total set of eigenvectors of Lirf 
It is clear from previous discussions that convergence of the iteration will depend on whether any of the nonzero 
eigenvalues has a magnitude greater than unity. 



SD 







" uC^) ' 




SD u^*^) 




" A.uC^) 


M 


T 








M uC'-) + T wC^-) 




AfewC^) 
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D. Fixed intermediate density and continuum models 



For additional insight into the the numerically relevant finite difference models of the previous section we now 
consider the continuum perturbation problem with 

Sk = -S^I^ . (3.32) 

It is convenient to view this in terms of the differential operator that is the inverse of the updating operator. As in 
Eq. ((2?2)) . we have 

1 d f 

and the eigenequation is 



L;^'{v)^{1/X)v. (3.34) 

The conditions on v{r) at r — > are as before, but now the other condition on the eigensolution is that v{rf) = 0. 
For the inner product 

/■'•/ 

Vi ■ V2 = / r'^'P{r)vi{r)v2{r) dr , (3.35) 
Jo 

this constitutes a Sturm-Liouville problem, and hence L^^ has a complete eigenbasis. More specifically, it has an 
eigenbasis that is complete in the mean (with respect to the above inner product) on the interval (0, r/). But the 
interval relevant to the stellar interior is (0, R), and there is no reason that the eigensolutions will be complete in any 
meaningful sense on (0, i?). 

The continuum eigenvectors for the interval (0,r/) correspond to the eigenvectors u'^'^^ of the (/ — 1) dimensional 
subspace in the discrete implementation of the problem. Just as the solutions to Eq. (|3.33p have no special significance 
for rf < r < R, the column eigenvectors v^*^^ of the discrete problem, in Eq. (|3.3ip . have no special significance for 
the bottom — / + 1 elements of the column. 

For rf < r < R, the continuum problem has similarities to the rf — fixed central density problem discussed in 
Sees. IIIIAI and HITB] In particular, 6{r - R) is an eigenfunction (or generalized function) with zero eigenvalue, and - 
as in Sees. IIII XI and UlI Bl - the generalized eigenvectors of the Jordan decomposition have analogs in the continuum. 
The physical picture of the zero eigenvector applies just as in Sees. IIIlXl and [Til Bl 

By specializing to the n — I case, we may once again benefit from a simple closed-form example. With the notation 
of Eq. (|2.6p . the operator (|3.23p becomes 

For this operator, 5{£^ — tt) is an eigenvector (in the vector space of distribution functions on R) with zero eigenvalue, 
and the solutions t;'^''"' to 

L,^(^W)=Afcz;W, (3.37) 

given by 

x;W(e)=r'sin(A;'/'0, A;^/' = fcTT/e/, fc=l,2,3---, (3.38) 

are eigenvectors. This suggests that, for an n = 1 polytrope, the iteration p.23p should converge for any choice of 
S,f < TT, but the convergence rate is maximized when the density is fixed at the center. We note, however, that the 
above eigenvectors are not complete on (0, tt) and that generalized eigenvectors can be constructed using the procedure 
of Eq. (IX^ . 
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IV. SUMMARY AND CONCLUSIONS 

We have investigated the properties of the updating operator for hnearized iteration of the two equations that 
govern Newtonian ncntron star structure, and have focused on spherically symmetric models. Wc have considered 
two constraints on the iteration: (i) the injection energy k is held fixed, and (ii) density is held fixed at a specified 
radius. 

In the case of fixcd-K iteration wc have found that both the finite dimensional problem (for finite difference discretiza- 
tion) and the continuum problem are self-adjoint and convergence is determined by the spectrum of its eigenvalues. 
For polytropic equations of state, numerical work had always led to divergence of iteration. We have shown, in 
fact, that there is a rigorous bound on the largest eigenvalue of the updating operator: it must be greater than the 
polytropic index n. Since all numerical experiments had been carried out with n > 1, this meant that there had to 
be an updating eigenvalue greater than unity, and hence that iteration would diverge. This did suggest, however, 
that for polytropic equations of state with unphysically small values of n, convergence might be possible. Numerical 
experiments showed that in fact this is the case, thereby verifying the applicability of the analysis. 

The updating operator for fixed central density was shown to have a very different nature than that for fixed-K. In 
the finite dimensional case corresponding to a finite difference representation of the equations, there is only a single 
eigenvector, with eigenvalue zero, and the updating operator is nilpotent. The continuum version of the fixed central 
density problem is not connected as directly to the finite difference problem as in the fixed-K case. We have shown, 
however, that for the continuum iteration converges. It is also possible to construct a sequence of functions that have 
some of the spirit of the generalized eigenbasis of the Jordan decomposition of the finite dimensional problem. 

For density fixed at some radius r/ other than the center, the updating operator, not surprisingly, has mixed 
properties. The space on which the updating operator acts can be separated into two sectors, one corresponding to 
radius ^ rf on which the updating operator acts more-or-less like the updating operator in the fixed-K case; and 
the other sector, for radius > r/ on which the updating operator acts more-or-less like the operator for fixed central 
density. 

Although only the finite difference analyses arc directly applicable to numerical iteration, the connection to the 
continuum is not only useful, but has a practical importance: it suggests that the properties of the updating operator 
are not idiosyncrasies of finite differences. It thus adds confidence that, for example, the use of a Gaussian method 
for the integral will converge or diverge just as the finite difference case would. 

The analyses presented here have been limited to spherical symmetry. But further work, mostly numerical, will 
be reported elsewhere that shows that many of the general conclusions reported here also apply to rotating stars 
and to binaries. In particular, the fixed-K updating operator is always self-adjoint, and generally has an eigenvalue 
greater than unity; and rotating stars with fixed rotation speed and fixed central density lead to a nilpotent updating 
operator. 

Although the main motivation for the work undertaken here has been a mathematical understanding of iteration 
properties, the results have potentially useful applications. In particular, the convergence of a nilpotent updating 
operator (like the fixed central density operator) is very different from that of a self-adjoint updating operator (like 
the fixed-K operator for a polytropic equation of state with a small polytropic index). The nilpotent operator will 
reach a solution at machine precision within a finite number of iterations, while convergent iteration in general may 
approach the correct solution gradually, and slowly. It may also be useful to understand the nature of the iteration 
when a test must be made whether or not the process is converging. For convergent iteration with a self-adjoint 
updating operator a simple measure of the difference in subsequent solutions can be used, employing the same metric 
for which the operator is self-adjoint. With an updating operator like that for fixed central density, convergence may 
give an early appearance of divergence. If the initial perturbation is a distribution concentrated near the outer edge 
of the stellar model, each iteration will move the perturbation inward, reducing it only after many steps. Too simple 
a test for convergence might misinterpret this iteration as nonconvergent. 
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